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Abstract. Using simulations and a simple mean-field theory, we investigate jamming 
transitions in a two-species lattice gas under non-equilibrium steady-state conditions. 
The two types of particles diffuse with different mobilities on a square lattice, subject 
to an excluded volume constraint and biased in opposite directions. Varying filling 
fraction, differential mobility, and drive, we map out the phase diagram, identifying 
first order and continuous transitions between a free-flowing disordered and a spatially 
inhomogeneous jammed phase. Ordered structures are observed to drift, with a 
characteristic velocity, in the direction of the more mobile species. 
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1. Introduction 

Since their introduction, driven lattice gases have attracted much attention as some of 
the simplest systems for exploring non-equilibrium steady states [TJ |2] . Motivated by 
the physics of fast ionic conductors [3], the earliest models consisted of a single species 
of "charged" particles, diffusing on a periodic lattice, subject to the effects of both a 
thermal bath and an external (uniform, DC) "electric" field. Driven by the field, the 
system settles into a state with non-trivial current. With no inter-particle interactions 
other than an excluded volume constraint, the stationary state distribution is trivial [3], 
and interesting behavior is displayed only in time- dependent quantities. However, once 
interactions are introduced, then a variety of unexpected phenomena arise even in the 
steady state, such as long-range correlations at all temperatures 0. 

A natural generalization is to study systems with two species, motivated by, e.g., 
fast ionic conductors with several mobile species [3J. Such systems have been used to 
model water droplets in microemulsions with distinct charges gel electrophoresis 
[3 |Hj, vacancy mediated diffusion in binary alloys and traffic flow [TU]. Focusing 
on systems with two species carrying equal but opposite "charge", we discovered 
remarkably complex steady states, even for particles subjected only to the excluded 
volume constraint [Hid]. ^ or simplicity, the early studies employed only square lattices 
of L x L, filled with equal numbers of "positive" and "negative" particles. The control 
parameters in this simple system are just E, the strength of the drive, and fh, the overall 
particle density (regardless of charge). By varying these parameters, we found that 
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this model displays a phase transition from a disordered state with homogeneous 
densities at small E or m, to an ordered state with inhomogeneous densities. Moreover, 
this line in the E-fh plane consists of a section with discontinuous transitions and another 
with continuous ones, which may be regarded as analogues of first and second order 
transitions in equilibrium systems. To understand these phenomena, it is sufficient to 
adopt a mean-field, continuum approach, based on hydrodynamic like equations for 
the two conserved densities. A linear stability analysis around the homogeneous state 
allows us to predict the existence of transitions. Fortunately, we were able to find 
analytic solutions for both, homogeneous and inhomogeneous steady states so that 
the presence of both types of transitions can be predicted. Subsequently, such models 
have been generalized to include rectangular lattices ^2], non- neutral ^3] or nearly 
filled ^Hj systems, and "charge exchange" dynamics ^H] in one- dimensional and 
quasi-one-dimensional fHJ EH] cases. A dazzling array of surprising phenomena was 
discovered. 

Here, we investigate a further generalization, namely, having species with different 
mobilities. Indeed, unless there are some underlying symmetry constraints, there is no 
reason to expect two different species of particles in typical physical systems to have 
identical mobilities. Using both Monte Carlo techniques and mean-field analyses, we find 
that, though the phase diagram appears to suffer little change, novel features arise in all 
regions of parameter space, even at not-so-subtle levels. In the next section, we present 
some details of how a differential mobility is implemented in simulations, followed by 
a discussion of our results. Section 3 is devoted to the continuum, mean-field theory 
and its predictions. We close with a summary of our findings and discuss some possible 
avenues for future investigations. 

2. Model Specification and Simulation Results 

Following the first studies of biased diffusion of two species, square lattices of size 
LxL = N are used. In most of our simulations, L = 40. We focus exclusively on neutral 
systems (i.e., N + = A_) at various densities fh = (N + + AL) /A. A configuration of 
the system is specified by the set {cr x>y }, where a X;V assumes the values 0, ±1 if the 
site (x, y) is empty or occupied by a positive (+) or negative (— ) particle. Clearly, 
?™ = J2x, y a x, y /N ■ The diffusive nature of the the particles is modeled by allowing them 
to jump only to nearest-neighbor empty sites. After a particle and one of its neighboring 
site are selected, a particle-hole pair is always exchanged, unless it results in a particle 
jumping against the external "electric" field. Chosen to point in the positive y axis, the 
field is parametrized by E = Ey, so that its sole effect is to lower the probability of 
particle jumps against E to e~ E . 

In contrast to previous studies, the two species are endowed with different mobilities, 
modeling say, trucks and cars on a road. Using the self-evident notation /i± for the 
mobilities, we arbitrarily choose /i + > /i_ and then focus on the ratio /i = To 
implement this ratio, we keep a list of the particles and their locations. In one Monte 
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Carlo step (MCS), fhN particles are chosen from the positive/negative list randomly 
with the frequency ratio of fx. In particular, by focusing on integer values of /x, we 
simply make fx attempts to move a positive particle for every attempt to move a negative 
one. By monitoring a few quantities, we found that there is little difference between 
sequential and random updating in species space, especially since the particle locations 
were randomly assigned at the beginning. After fhN attempts, we increment the time 
by 1 MCS. Below, we will see that it is often convenient to use the "renormalized" MCS 
(1 RMCS being defined as fx + 1 MCS), during which every slow particle would have 
had, on the average, one chance to move. Typically, 20K MCS or more are discarded 
before measurements are taken. In this study, we probe fx — 1,3, 10, 30, 100, and 300, 
though most of the data are collected for fx — 30, being a compromise between exploring 
large differentials and dealing with finite computation times. 



2.1. Continuous and discontinuous phase transitions 

Similar to the fx — 1 system, ours displays a transition from a disordered homogeneous to 
an ordered inhomogeneous state, the latter showing a marcoscopic cluster with opposing 
particles locked in a "jam." To differentiate these phases, we monitor the first few mass 
structure factors, defined as 

S n = ^|?7Jn| ^ , 

where 

fh n = L~ 2 J2 a l, y ex P [2mny/L] 

x,y 

is the Fourier tranform of the mass density profile across y. Here, (•) denotes the average 
over 800 configurations, separated from one another by 100 MCS in a typical run of 
80K MCS. We also measure the fluctuations in the S"s by monitoring the variances: 
^ \m n \ 2 — S n ^. Both continuous and discontinuous transitions are observed, as 

illustrated by a phase diagram in the E-fh plane for the fx — 30 case (Fig. 1). The 
insets illustrate the nature of the transitions, in two typical parameter domains. 

To provide further detail on how these lines are determined, we expand the insets 
of Fig. 1 into Fig. 2a,b. Here, we plot the variations in S 1 as the system is swept back 
and forth across the transitions. For the continuous case (Fig. 2a), we show the result 
from a sweep in £ 6 [0, 0.5] while fh is held at 0.6. This direction of traversing the 
transition line is chosen because this line is almost parallel to the m-axis in this region. 
Note that there are two sets of data points (x/— ), representing the up/down sweeps 
from a single long run. This run is started with E = in a random configuration at 
just over half filling (fh = 0.6). Then, the first 20K MCS are discarded, data gathered 
for the next 80K MCS, E raised by 0.02, 20K discarded and so on, until E = 1.0. 
Thereafter, the process is reversed, decreasing E until E — 0. The figure shows clearly 
that the transition is continuous with no sign of hysteresis. At the crude level of our 
investigations, the critical field is then identified through either the largest dSi/dE 
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Figure 1. Phase diagram for mobility ratio of 30, with the system being homogeneous 
and "jammed" in the lower-left and upper-right regions, respectively. The solid line 
marks continuous transitions. For discontinuous transitions, a pair of symbols (+/— ) 
are used to denote the jumps in the order parameter when the filling fraction (fh) 
is increased/decreased. The absence/presence of hysteresis loops is illustrated by 
the insets, each containing data for sweeps in both directions (cf. Fig. 2 for further 
explanations). 



or the peak in the variance of S\. The transition line is assembled by repeating such 
runs for fh G [0.5, 0.95] in steps of 0.05. For the discontinuous cases, which are more 
easily probed by sweeping in fh, we observe the typical hysteresis loops. The inset in 
Fig. 1 and Fig. 2b show the case of E — 2. Again, the data represent a single long 
run, starting with m = 0.1 in a random configuration, discarding the first 20K MCS, 
gathering data for the next 80K MCS as above, raising fh by 0.01 and so on. These 
data points are represented by crosses (x). After a run at fh — 0.4, we continue this 
process with decreasing m's until fh = 0.1. This set of data is shown as minus signs 
(— ). A hysteresis loop is clearly displayed, and the densities at the jumps are recorded 
to construct the points plotted in the main figure of Fig. 1. For comparison with the 
previous studies (i.e., the \x = 1 case), we also show the results of such runs as solid lines 
in Fig. 2a,b. As we see, there is little difference between the \x = 30 and /j, — 1 results, 
leading us to conclude that, within the accuracy and statistics of our study, there is no 
significant dependence of the phase diagram on \x. 

Before turning to the properties of the inhomogeneous state, let us remark that, of 
course, there are non-trivial effects due to \i > 1, even in the disordered phase. Since 
one species is more mobile, there must be a difference between the average velocities, 
which in turn leads to a non-zero mass current even though iV + = AL. To verify this 
expectation, we measure the two currents, J±, by simply keeping track of the number 
of +/— particles which move up/down. Since iV + = AL in our samples, the ratio of 
the average velocities is just v + /v- = J+/J-. With limited observations, our naive 
expectation: 

J+ / J- = n 
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Figure 2. Plots of Si as the system is swept across the transitions, (a) E is increased 
(x) and then decreased (— ), with fh — 0.6 and fi = 30. (b) m is increased (x) and 
then decreased (— ), with E = 2 and [i = 30. For comparison, the /i = 1 cases are 
shown as solid lines in both panels. With our statistics and accuracy, the effects of 
mobility differential are relatively minor. 

turns out to be satisfied quite well. As the data in Table 1 show, we have investigated 
only a few sample points, in both "arms" of the disordered region (small E moderate 
m, and vice versa). 
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Table 1. Currents and their ratios for various regions in the homogeneous phase. 

Naturally, we expect other interesting properties associated with the disordered 
state; yet, these are generally quite subtle. For example, even in the \i = 1 case, there 
are singular structure factors and long-range correlations |2Hj • For more prominent and 
fascinating features, we turn to the inhomogeneous state. 

2.2. Drifting structures in the ordered state 

Here, the density of particles is high enough so that they "jam" into a macroscopic 
cluster. Consequently, the translational symmetry (in y) is spontaneously broken. In 
neutral systems with \i — 1, the cluster performs an unbiased random walk, and, due 
to the symmetry under charge "conjugation" (+ <^ — ), its charge profile is purely 
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Figure 3. Drift speed as a function of /i. The inset shows the "center of mass" of the 
macroscopic cluster, as a function of time, for three values of [i. 



antisymmetric. On the other hand, for charged systems (i.e., iV + 7^ AL) this cluster 
drifts " backwards" (opposite to the average motion of the majority species)! Though 
counter-intuitive at first glance, this behavior can be understood [Tl|. Simulating > 1 
systems here, we find the cluster to drift "forwards," i.e., in the direction of the more 
mobile species. Though this behavior may seem less surprising, it is easy to advance a 
(wrong) heuristic argument to arrive at the opposite conclusion! 

To characterize the drifting cluster more systematically, we carry out the following 
analysis. Since the runs are necessarily long, we focus only on neutral lattices with 
fh = 0.36, driven with E = 1. These parameters put us well in the ordered phase. To 
study the /i-dependence of the drift speed, we perform simulations with /1 e [2, 32] in 
steps of 2. For comparison with earlier studies, we also collect data for jj, — 1. Using a 
random configuration at the start, the usual first 20K MCS are discarded. To be sure 
that a macroscopic cluster is present, we first test the magnitude \rhi\ and proceed only 
when it is larger than 0.25 (c.f., \rhi\ = 1 /[L sin(7r/L)] ~ 0.32 for a fully ordered state). 
Assuming the center of mass (CM) of the cluster is well characterized by the CM of the 
entire configuration, we measure the latter at intervals of 500 RMCS. Here, the CM of 
a particular configuration is defined simply by the phase in rhi = \rhi\e 1,9 . Also, we find 
it more convenient to use the RMCS as a unit, since the time scale of the system as a 
whole is really set by the slow-moving particles. In the inset of Fig. 3, we show how 9 
depends on time for three /z's. Clearly, the dominant behavior is a steady drift. Fitting 
these data to vet + const, we extract the drift speed vq. Translating vg into a velocity, v, 
in units of lattice spacing per 1000 RMCS, we plot v (n) in Fig. 3. Until we have some 
understanding of this quantity, we refrain from fitting it to a particular form. Instead, 
we just note that, interestingly, v seems to saturate at unity. Finally, we have computed 
the standard deviations from the linear fits naively and found that these are essentially 
independent of fi. Of course, since we believe the CMs to be performing (biased) random 
walks, we should expect these deviations to grow as t x l 2 , from which diffusion coefficients 
can be extracted. Postponing such a detailed study to a later publication, what we may 
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Figure 4. Mass (a) and charge (b) profiles in the ordered state. The dashed line and 
the open circles correspond to the y, = 1 and y = 30 cases, respectively. Note that the 
charge profile for the latter is not purely antisymmetric about the CM (L/2 = 20). 



infer from our "naive" computation is just that clusters associated with different /x's 
execute random walks with the same diffusion constant but different biases. 

Next, we consider the mass and charge profiles in the jammed state. For each 
configuration we measure both 6 and 

m (y) = S a ly and i(y) = Yl a ^v ■ 

y y 
Then we displace the y-dependent quantities according to y — > u = y — 9L/2tt + L/2, so 
that the CM is located at L/2 = 20 for convenience. In this co- moving frame, we can 
now average over the run to produce the profiles: (m (u)) and {q {u)). Plotted in Fig. 4 
are these profiles for the parameter set (E,m,fi) = (1,0.36,30), and, for comparison 
with previous results, those for the ji = 1 case as well. Even though the change 
in the mass profile appears almost imperceptible, there are actually some significant 
differences. For example, though the densities within (outside) the cluster differ by only 
about +0.003 (—0.003), this discrepancy translates into a factor of two for the densities 
far from the "jam." Until we have much better statistics, it is unclear if we can attach 
much meaning to these subtle differences. The changes in the charge profile are more 
discernable. We may interpret these as an enhancement of the more mobile species' 
ability to penetrate into the block of the less mobile one. 



3. Mean-Field Theory and Comparisons 



To develop a first-step understanding of the large scale phenomena shown above, we 
exploit the same mean field approach which proved quite successful in the fi = 1 case 
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P"H IT3] . In this approach, hydro dynamic like equations of motion are postulated for the 
densities of the two species, p± (x, t) = (a 2 (x, t) ± a (x, t)) /2 : 

where J± are the associated (density dependent) currents [TT]. Here, we also treat 
the spatial co-ordinates as continuous variables, in general d dimensions: x = 
(xi, ...,Xd = y). To model differential mobility, we simply multiply the currents by p±. 
Absorbing the average mobility into the time scale by defining 

M± EE (1 ± 8) , 

the continuum equations read 



dp± 
dt 

where 



i±«j)v- </>Vp±-p±v<pTP±<f£y , (i) 



<t> = i - p+ - p- 

is the hole density and S represents the coarse-grained electric field {£ = 2 tanh (£7/2) 
at the most naive level). Of course, these equations can be re-expressed in terms of 
and the charge density 

V> = P+~ P- 
resulting in a form 



(2) 



d_( \ [ 1 -5]( V-[V0 + ^y] \ 

dt\ 4> J [-5 1 J V V ■ \<pVi) - i)V<p - 0(1 - 4>)Sy] J 

which easily reduces to the one studied previously [TTJ [TH] . 

To find solutions to these equations, we must supply boundary conditions and 
constraints. For a hypercubic system of volume L d , they are 

P± (x, t) = p± (x + Lxi, t) , i = 1, 2, d 

and 

dx _ f dx 

—j<p{x,t) = l-m, J j^if){x,t) = 0. 

We will discuss mainly steady states, i.e., density profiles which are time independent 
in a co-moving frame. From Monte Carlo results as well as previous experience with 
this model, we expect these profile to obey translational invariance in the transverse 
directions in the whole parameter space. On the other hand, spontaneous breaking of 
translational invariance in y is manifest for the jammed state. Due to the drift, we will 
denote these profiles by (ft (u) and ip (u), with u — y — vt. 
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3.1. Homogeneous states and linear stability analysis 

Let us consider first the disordered phase, where the steady-state densities are 
homogeneous in both space and time, i.e., 

4> (x, t) — 1 — m, if; (x, t) = 0. 

Note that, unlike the \l — 1 case, the mass current is no longer zero. The (average), 
denoted by Cy, is given through 

C = y • (J+ + J_) = — m)m/2 — — m)m/2 = 5rh(l — rh)S. 

With our normalization of time scales, the charge current is unchanged: 

J = y ■ (^J + — J_j = m(l — m)E. 

Some care is needed when direct comparisons of absolute currents with Monte Carlo 
data are made, since MCS and RMCS differ by a factor of fi + 1. One test of this theory 
which is independent of such details is the ratio of current magnitudes 

J+/J- = /i. 

As noted above, this prediction agrees well with our limited observations (apart from 
the first set, in which both currents are very small and noisy). Of course, our crude 
data should be regarded as merely the first step in a more refined study. That they are 
consistent with J+/J- = /U just confirms that this naive way of introducing mobility 
differentials is a reasonable starting point. Noting that the agreement is somewhat worse 
for high densities suggests the presence of subtle particle-particle correlations, which can 
only be uncovered by better statistics. 

Beyond these trivial results, a linear stability analysis around this state leads us to 
the matrix 

— k 2 ik y £(l — rh) 

ik y S(2rh-l) -k 2 {l-m) _ 

which is purposely written in a form to display the effects of 5 > 0. In contrast to the 
6 = case, the eigenvalues (Ai^) of this M are generically complex. Thus, we should 
examine their real parts to identify the stability limit. However, since M is real, it 
can be expressed as |Ai| 2 Re A2/ Re Ai (or 1 <^ 2), so that detM = is still valid for 
identifying the stability limit. Since det M is modified by a simple factor of (1 — 5 2 ), we 
conclude that differential mobility has no effect on this line! In other words, the most 
unstable perturbation is still given by the lowest longitudinal mode: 

k y = 27r/L 

and the instability occurs at [TT] 

(£L/2itf (2m - 1) = 1. 

Of course we cannot expect good quantitative comparisons with the phase boundaries 
from Monte Carlo simulations. For example, in the Ising model, mean-field analysis 
overestimates the critical temperature by nearly a factor of 2. Nevertheless, we are 
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encouraged by one aspect of this prediction, namely, that differential mobility has no 
effect on the onset of the instability. This seems to be consistent with the observed 
phase boundaries being quite insensitive to //. 

For completeness, we give the explicit expressions for the eigenvalues: 

A 12 = -k 2 (l - p) - 5ik y £p ± VR (3) 



where 



R 



p 2 + 5 2 (1 - m)\ - (k y S) 2 [(1 - 5 2 )(2m - 1)(1 - m) + 5 2 p 2 

+ 25ik 2 k y £p(l - p) (4) 

and p = m/2. Their complex nature reflects drifts associated with the eigen- 
perturbations. Given that the species are mobile in different ways, such drifts are 
to be expected. In principle, it is possible to measure the decay and drift of small 
perturbations in simulations so that comparisons with these predictions can be made. 
In practice however, such measurements would require better statistics and longer runs 
than those in this study. 

3.2. Density profiles in inhomogeneous states 

Turning to the more interesting ordered states, the most immediate question concerns 
the steady-state profiles. As observed in simulations and expected from symmetry 
considerations, these profiles are homogeneous in x but depend on y. In analogy to 
the /x = 1 case, we seek non-trivial functions 0* (y) and ip* (y) which satisfy Eqn. (j2J) 
in a steady state. However, it is clear from Eqn. (0) that if non-trivial functions lead 
to zeroes on the right hand side with /i = 1, then they also lead to zeroes for any p). 
The conclusion is that our mean-field theory admits inhomogeneous profiles which are 
(i) stationary (i.e., do not drift) and (ii) identical to those for p, = 1. Of course, we may 
argue that this obvious discrepancy, when comparing with data, is a "small" effect in 
absolute terms. However, our equations are non-linear and the existence of additional, 
drifting solutions cannot be ruled out. Unfortunately, we are unable so far, either 
to find drifting solutions (analytically or numerically), or to prove that our equations 
admit no such solutions. Nevertheless, for completeness, let us present the analysis 
which simplifies the problem to a single second order, non-linear, ordinary differential 
equation. 

Following the techniques in ^3], we assume the forms 
/ <f>*(x,t) 

where 

u = y — v t, 

with v being the drift velocity of the macroscopic cluster. Inserting these into Eqn. (j2J, 
we obtain 

[(d(f)*/du) + <p*i)*E] 
[(f)* (d^*/du) - i)* {df /dv) - <f(l - (jf)S\ 




\ d 


1 -6' 


1 du 


-8 1 
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which has the form of a vanishing total derivative. So, the first integrals are just 
constants which are identified with the two currents. As presented above, we denote the 
mass current C and the charge-current by J. Recognizing that the hole current must 
be — C, we obtain 



(5) 




' [{d<p*/du) + </)*iIj*£} 

v [0* (d^*/du) - i)* (d(j)*/du) - 0* 

To simplify these equations, we change variables, as in the 5 = case, to 

X=l/(f>* and u = ip*/(p*. 
Then, we invert the matrix so that the first derivatives are decoupled: 



1 - 0* )£] 



(—dx/du + ujS 
duo/du — (x ~ 1)£ 



-x 



1 5 
5 1 




1-5 2 

Being first order differential equations, their full solutions will require two more unknown 
integration constants. Together with C and J, there are four unknowns to be fixed, by 
the four constraint equations: 

X(0)= X (L) and u (0) = cu (L) 

from periodic boundary conditions, as well as 

f -du= (<f>* = L{l-m) and f - = f ^* = 
jo x J J X J 

from the conservation laws. Further simplifications occur when we define C, J, v with 

appropriate factors of S, 5, and 1 — 5 2 : 

C 

= ssc 

= £J 



1-5 2 
J 

V 



1-5 2 

so as to exploit the scaling of u to the variable u£, and the symmetry of the system 
under 5 — 5 (with C, J, v even in 5). Finally, denoting d/d{u£) by prime ('), we arrive 
at 



X ' = uj + 5(J-C) x 2 + Sv X (l + Su) 
uj'=( X -1)-(J-S 2 C) X 2 -Sv X (5 



UJ 



(6) 
(7) 



One advantage of this form lies in that, in the limit 8 — > 0, it easily reduces to the 
previous case: x" — (x ~ 1) — Jx 2 - The other is the clear separation of x, u; into 
even/odd functions of 5. Note that, under this "parity" operation, u£ is odd. 

As in the 5 = case, the first equation can be trivially solved for ui and the resultant 
inserted into the second, leading us to a single (though quite complicated) equation for 
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X- As pointed out above, there exists at least one solution for any 5, corresponding to 
a stationary (v = 0) state and involving C = J (i.e., J = 5C). 

Although we have not found non-trivial drifting solutions analytically, we are able 
to test the consistency of these equations against the measured profiles and v. For 
example, we can begin with Eqn.fjHJ) and integrate over u to get relationships between 
v and averages of the profiles. Indeed, there are infinitely many similar such relations, 
obtained from first multiplying this equation by functions of \ (or <p*) before integration. 
To minimize possible large errors, we choose to consider two equations so that the 
averages involve no higher power of x than unity. Specifically, we divide Eqn.fjHJ) by \ 
and x 2 before integration. The results are 

= (ijjcf)) +8(j-C)+8v(l-fh) 

= 5(J-C) ( x ) + 5v(l + 5{u)) 

where the average is defined via 




Eliminating the unknown currents, we arrive at 

y/s m (x) 

1 _ 6 2 j _ (0) {x) + s {rx) 

To compare with data, we retrace the rescaling of time and re-introduce the two 
mobilities, so that 

2//+/X- (X) vtf)l 

v = . 

(fi+ + //_) [1 - (1 - m) (x)} + (//+ - y.-) (i(f*x) 

Setting /x + = 30 and //_ = 1 would correspond to using RMCS as a unit of time. Using 

the naive £ = 2tanh(£'/2) with E = 1 and computing the appropriate averages from 

the measured profiles (Fig. 4), this relation predicts v ~ 0.004. Being both of the right 

sign and order of magnitude, this is an encouraging result. 

There is an alternative approach, based on a discrete version of a simple hopping 

model. Without delving into the details, we only state that it can be obtained as an 

approximation from the full master equation. By considering (a {x,t)) and ignoring all 

correlations, the result is a discrete equation for the average densities p± (x, t). Focusing 

only on the y co-ordinate (i.e., averaging the densities over x), our starting point is 

P± (y, t + l)-p± (y, t) = i Q { J in - J out } , (8) 

where the jumps into/out-of site y are given by 

Jin = P± (y =F 1, t) (j) (y, t) + e~ E p± (y±l,t)(f> (y, t) 
Jout = p± (y, t)(p(y± 1, t) + e~ E p± (y, t) (y =F 1, t) . 

Here the time unit is a RMCS, while the factor 1/4 accounts for half of jumps being 
transverse, so that each term in the J's represents jumps with 1/4 probability. Of 
course, periodic boundary conditions are imposed: p±(y,t) = p± (y + L,t). For an 
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inhomogeneous drifting steady state, we seek functions of the form p± (y, t) 
so that we would write naively 



Pi w 



P*± (u-v)- p± (u) 



+ e 



P± («Tl) 
p* ± (u±l) 



u 



P± w 
u) - p* ± (it) 



u±l) 
«=Fl)]} 



Expecting v to be of the order of 10~ 3 , the first term on the left must be modified if 
we wish to apply this equation to data analysis. We believe a good approximation is to 
interpolate the left hand side: 

P*± ( u ~ V )-P*± 0) - v [p± ( u ~ !) - P± ( u ) 
Rearranging the result, we have 

4v [pi (u - 1) - pi («)] = (A {e~ E p* ± (u ± 1) 0* (u) - p* ± (u) <P* [u ± 1) 

+ P * ± (ut 1) 0* (u) - e~ E p* ± (it) 0* (u T 1)} 

which can be "integrated" once, as usual. The constants are just the steady state 
currents (J± = y ■ J±), so that our basic equations read 

(9) 
(10) 



4J; = p [p* + (u) 0* (u + 1) - e~ E p\ (u + 1) 0* (u)\ - Avp\ (u) 
AJ*_ = \e~ E p*_ (u) 0* (u + 1) - pi (w + 1) 0* (uj\ - Avp*_ (u) 



In principle, these are recursive relations for the profiles and, imposing constraints 
(J2 p± = N±), the unknown profiles and constants ( and v ) can be found. In practice, 
this method is not simple. Here, let us simply test these relations against the profiles 
from the data. By regarding e~ E also as an unknown, we can write four linear equations 
for these constants and see how well the agreement is. For simplicity, in a manner similar 
to the continuum study above, we choose to sum Eqns. ()9I10|) for two equations. For 
the other pair, we first divide Eqns. ()9|1()|1 by p* ± {u) respectively and then perform the 
sum. The results are 

0.99 x 10~ 2 
-0.65 x 10~ 3 
0.63 x 10" 3 
0.39 

Again, we are encouraged by how well this crude scheme functions, since the measured 
v is about 0.9 x 10~ 3 and e _1 ~ 0.37. As for the currents, their small values are typical 
of jammed states. Though our data for them are too noisy for a meaningful comparison, 
we are quite satisfied by the relative magnitudes and signs. 



J* 




J* 




V 









4. Summary and Outlook 



We have reported simulation results and mean-field arguments for jamming transitions in 
a three-state lattice gas under non-equilibrium conditions. In our simple model, positive 



Effects of differential mobility on biased diffusion of two species 



14 



and negative particles diffuse on a lattice, subject to an excluded volume constraint and 
an external drive which biases their motion in opposite directions. Generalizing earlier 
studies of this model^Tj, we focus here on the effect of differential mobility, allowing one 
species to be more mobile than the other by a factor of //. With equal numbers of the 
two species in the system, our key findings are as follows: The previously observed phase 
transition persists for the range \i G [1,300], i.e., a free-flowing state, for low densities 
ifh) and drive (E), giving way to a jammed phase at higher rh and E. Moreover, 
there is little quantitative change to the fh-E phase diagram. Crossing the transition 
line at low rh and high E, we observe hysteresis loops in the density structure factor, 
consistent with first order transitions. Across the high rh and low E portion of the line, 
the transitions are continuous. While the phase diagram is essentially unaffected by 
\x > 1, both disordered and ordered phases exhibit a systematic drift, i.e., non- vanishing 
mass current. Focusing on ordered structures, we find that the cluster drifts in the 
direction of the more mobile species. According to the data, hole and charge density 
profiles depend in subtle ways on /i, and their drift velocity appears to saturate as 
\x increases. To shed more light on these findings, we present a mean-field theory, in 
the form of coupled, nonlinear partial differential equations of motion for the hole and 
charge densities. Homogeneous solutions, corresponding to disordered phases, follow 
trivially from the conservation laws for the two particle numbers. A linear stability 
analysis demonstrates the presence of an instability, as the overall density increases. 
Unfortunately, we were not able to find inhomogeneous drifting solutions, so that we 
can only offer several consistency checks between equations and data, in order to build 
confidence in our theoretical description. 

Many questions remain open for further study. First of all, a better understanding of 
the drift velocity would be desirable. Are there inhomogeneous solutions to our mean- 
field equations with non-trivial drift? If a proof of their absence can be established, 
then we should include noise terms (to promote the mean-field equations to Langevin 
equations) and study the effects of fluctuations and correlations. Perhaps these will 
renormalize the "bare" co-efficients in Eqn. (JI} so as to admit drifting structures. 
Such analysis, along with further simulations, should also settle the question whether 
the velocity continues to grow or approaches a saturation value as /i increases. In 
addition, the noise terms will allow us to investigate particle correlations and structure 
factors. As shown in preceding studies |20j . these can be extremely interesting, even in 
the disordered phase. Going beyond equal-time structure factors, a study of dynamic 
correlations would provide considerable insight into currents and drifts. Finally, we note 
that, at the first order line, we should find coexisting phases. Since the mass current is 
significant in the disordered region but quite small in the jammed region, these states 
will undoubtedly display a rich variety of behavior. 

Nearly all the simulations performed here are based on a 40 x 40 lattice. For fixed 
E, the observed transitions cannot possibly persist as the longitudinal size (L y ) goes to 
infinity. At the mean-field level, the inhomogeneous solutions effectively depend only 
on the product £L y . Thus, performing simulations with a range of E and L y would be 
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useful for establishing the presence of u EL y scaling" or proving its absence. In general, 
explorations of finite size effects are clearly crucial before reliable conclusions on the 
nature of the phases and transitions can be drawn. 

A natural expansion of parameter space is to allow for unequal numbers of the 
species. Such systems exhibit drifting structures even in the absence of differential 
mobility[14 . As a consequence, one might wonder whether it is possible to adjust 
these particle numbers and \i such that the jam becomes stationary. Viewed as traffic 
problems, involving cars and trucks at different densities and mobilities, the answers to 
these questions may shed some light on the essence of traffic jams. 
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